home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / qr.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  13.2 KB  |  545 lines
  1. /* linalg/qr.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <string.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_vector.h>
  27. #include <gsl/gsl_matrix.h>
  28. #include <gsl/gsl_blas.h>
  29.  
  30. #include <gsl/gsl_linalg.h>
  31.  
  32. #define REAL double
  33.  
  34. #include "givens.c"
  35. #include "apply_givens.c"
  36.  
  37. /* Factorise a general M x N matrix A into
  38.  *  
  39.  *   A = Q R
  40.  *
  41.  * where Q is orthogonal (M x M) and R is upper triangular (M x N).
  42.  *
  43.  * Q is stored as a packed set of Householder transformations in the
  44.  * strict lower triangular part of the input matrix.
  45.  *
  46.  * R is stored in the diagonal and upper triangle of the input matrix.
  47.  *
  48.  * The full matrix for Q can be obtained as the product
  49.  *
  50.  *       Q = Q_k .. Q_2 Q_1
  51.  *
  52.  * where k = MIN(M,N) and
  53.  *
  54.  *       Q_i = (I - tau_i * v_i * v_i')
  55.  *
  56.  * and where v_i is a Householder vector
  57.  *
  58.  *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
  59.  *
  60.  * This storage scheme is the same as in LAPACK.  */
  61.  
  62. int
  63. gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
  64. {
  65.   const size_t M = A->size1;
  66.   const size_t N = A->size2;
  67.  
  68.   if (tau->size != GSL_MIN (M, N))
  69.     {
  70.       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  71.     }
  72.   else
  73.     {
  74.       size_t i;
  75.  
  76.       for (i = 0; i < GSL_MIN (M, N); i++)
  77.     {
  78.       /* Compute the Householder transformation to reduce the j-th
  79.          column of the matrix to a multiple of the j-th unit vector */
  80.  
  81.       gsl_vector_view c_full = gsl_matrix_column (A, i);
  82.           gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
  83.  
  84.       double tau_i = gsl_linalg_householder_transform (&(c.vector));
  85.  
  86.       gsl_vector_set (tau, i, tau_i);
  87.  
  88.       /* Apply the transformation to the remaining columns and
  89.          update the norms */
  90.  
  91.       if (i + 1 < N)
  92.         {
  93.           gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
  94.           gsl_linalg_householder_hm (tau_i, &(c.vector), &(m.matrix));
  95.         }
  96.     }
  97.  
  98.       return GSL_SUCCESS;
  99.     }
  100. }
  101.  
  102. /* Solves the system A x = b using the QR factorisation,
  103.  
  104.  *  R x = Q^T b
  105.  *
  106.  * to obtain x. Based on SLATEC code. 
  107.  */
  108.  
  109. int
  110. gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
  111. {
  112.   if (QR->size1 != QR->size2)
  113.     {
  114.       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
  115.     }
  116.   else if (QR->size1 != b->size)
  117.     {
  118.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  119.     }
  120.   else if (QR->size2 != x->size)
  121.     {
  122.       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
  123.     }
  124.   else
  125.     {
  126.       /* Copy x <- b */
  127.  
  128.       gsl_vector_memcpy (x, b);
  129.  
  130.       /* Solve for x */
  131.  
  132.       gsl_linalg_QR_svx (QR, tau, x);
  133.  
  134.       return GSL_SUCCESS;
  135.     }
  136. }
  137.  
  138. /* Solves the system A x = b in place using the QR factorisation,
  139.  
  140.  *  R x = Q^T b
  141.  *
  142.  * to obtain x. Based on SLATEC code. 
  143.  */
  144.  
  145. int
  146. gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
  147. {
  148.  
  149.   if (QR->size1 != QR->size2)
  150.     {
  151.       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
  152.     }
  153.   else if (QR->size1 != x->size)
  154.     {
  155.       GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
  156.     }
  157.   else
  158.     {
  159.       /* compute rhs = Q^T b */
  160.  
  161.       gsl_linalg_QR_QTvec (QR, tau, x);
  162.  
  163.       /* Solve R x = rhs, storing x in-place */
  164.  
  165.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
  166.  
  167.       return GSL_SUCCESS;
  168.     }
  169. }
  170.  
  171.  
  172. /* Find the least squares solution to the overdetermined system 
  173.  *
  174.  *   A x = b 
  175.  *  
  176.  * for M >= N using the QR factorization A = Q R. 
  177.  */
  178.  
  179. int
  180. gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
  181. {
  182.   const size_t M = QR->size1;
  183.   const size_t N = QR->size2;
  184.  
  185.   if (M < N)
  186.     {
  187.       GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
  188.     }
  189.   else if (M != b->size)
  190.     {
  191.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  192.     }
  193.   else if (N != x->size)
  194.     {
  195.       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
  196.     }
  197.   else if (M != residual->size)
  198.     {
  199.       GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
  200.     }
  201.   else
  202.     {
  203.       gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
  204.       gsl_vector_view c = gsl_vector_subvector(residual, 0, N);
  205.  
  206.       gsl_vector_memcpy(residual, b);
  207.  
  208.       /* compute rhs = Q^T b */
  209.  
  210.       gsl_linalg_QR_QTvec (QR, tau, residual);
  211.  
  212.       /* Solve R x = rhs */
  213.  
  214.       gsl_vector_memcpy(x, &(c.vector));
  215.  
  216.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);
  217.  
  218.       /* Compute residual = b - A x = Q (Q^T b - R x) */
  219.       
  220.       gsl_vector_set_zero(&(c.vector));
  221.  
  222.       gsl_linalg_QR_Qvec(QR, tau, residual);
  223.  
  224.       return GSL_SUCCESS;
  225.     }
  226. }
  227.  
  228.  
  229. int
  230. gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
  231. {
  232.   if (QR->size1 != QR->size2)
  233.     {
  234.       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
  235.     }
  236.   else if (QR->size1 != b->size)
  237.     {
  238.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  239.     }
  240.   else if (QR->size2 != x->size)
  241.     {
  242.       GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
  243.     }
  244.   else
  245.     {
  246.       /* Copy x <- b */
  247.  
  248.       gsl_vector_memcpy (x, b);
  249.  
  250.       /* Solve R x = b, storing x in-place */
  251.  
  252.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
  253.  
  254.       return GSL_SUCCESS;
  255.     }
  256. }
  257.  
  258.  
  259. int
  260. gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
  261. {
  262.   if (QR->size1 != QR->size2)
  263.     {
  264.       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
  265.     }
  266.   else if (QR->size1 != x->size)
  267.     {
  268.       GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
  269.     }
  270.   else
  271.     {
  272.       /* Solve R x = b, storing x in-place */
  273.  
  274.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
  275.  
  276.       return GSL_SUCCESS;
  277.     }
  278. }
  279.  
  280. int
  281. gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
  282. {
  283.   if (R->size1 != R->size2)
  284.     {
  285.       GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
  286.     }
  287.   else if (R->size1 != b->size)
  288.     {
  289.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  290.     }
  291.   else if (R->size2 != x->size)
  292.     {
  293.       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
  294.     }
  295.   else
  296.     {
  297.       /* Copy x <- b */
  298.  
  299.       gsl_vector_memcpy (x, b);
  300.  
  301.       /* Solve R x = b, storing x inplace in b */
  302.  
  303.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
  304.  
  305.       return GSL_SUCCESS;
  306.     }
  307. }
  308.  
  309.  
  310. /* Form the product Q^T v  from a QR factorized matrix 
  311.  */
  312.  
  313. int
  314. gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
  315. {
  316.   const size_t M = QR->size1;
  317.   const size_t N = QR->size2;
  318.  
  319.   if (tau->size != GSL_MIN (M, N))
  320.     {
  321.       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  322.     }
  323.   else if (v->size != M)
  324.     {
  325.       GSL_ERROR ("vector size must be N", GSL_EBADLEN);
  326.     }
  327.   else
  328.     {
  329.       size_t i;
  330.  
  331.       /* compute Q^T v */
  332.  
  333.       for (i = 0; i < GSL_MIN (M, N); i++)
  334.     {
  335.       gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
  336.           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
  337.       gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
  338.       double ti = gsl_vector_get (tau, i);
  339.       gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
  340.     }
  341.       return GSL_SUCCESS;
  342.     }
  343. }
  344.  
  345.  
  346. int
  347. gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
  348. {
  349.   const size_t M = QR->size1;
  350.   const size_t N = QR->size2;
  351.  
  352.   if (tau->size != GSL_MIN (M, N))
  353.     {
  354.       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  355.     }
  356.   else if (v->size != M)
  357.     {
  358.       GSL_ERROR ("vector size must be N", GSL_EBADLEN);
  359.     }
  360.   else
  361.     {
  362.       size_t i;
  363.  
  364.       /* compute Q^T v */
  365.  
  366.       for (i = GSL_MIN (M, N); i > 0 && i--;)
  367.     {
  368.       gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
  369.           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), 
  370.                                                                 i, M - i);
  371.       gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
  372.       double ti = gsl_vector_get (tau, i);
  373.       gsl_linalg_householder_hv (ti, &h.vector, &w.vector);
  374.     }
  375.       return GSL_SUCCESS;
  376.     }
  377. }
  378.  
  379.  
  380. /*  Form the orthogonal matrix Q from the packed QR matrix */
  381.  
  382. int
  383. gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
  384. {
  385.   const size_t M = QR->size1;
  386.   const size_t N = QR->size2;
  387.  
  388.   if (Q->size1 != M || Q->size2 != M)
  389.     {
  390.       GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
  391.     }
  392.   else if (R->size1 != M || R->size2 != N)
  393.     {
  394.       GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
  395.     }
  396.   else if (tau->size != GSL_MIN (M, N))
  397.     {
  398.       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  399.     }
  400.   else
  401.     {
  402.       size_t i, j;
  403.  
  404.       /* Initialize Q to the identity */
  405.  
  406.       gsl_matrix_set_identity (Q);
  407.  
  408.       for (i = GSL_MIN (M, N); i > 0 && i--;)
  409.     {
  410.           gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
  411.           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
  412.                                                                 i, M - i);
  413.           gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
  414.           double ti = gsl_vector_get (tau, i);
  415.           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
  416.     }
  417.  
  418.       /*  Form the right triangular matrix R from a packed QR matrix */
  419.  
  420.       for (i = 0; i < M; i++)
  421.     {
  422.       for (j = 0; j < i && j < N; j++)
  423.         gsl_matrix_set (R, i, j, 0.0);
  424.  
  425.       for (j = i; j < N; j++)
  426.         gsl_matrix_set (R, i, j, gsl_matrix_get (QR, i, j));
  427.     }
  428.  
  429.       return GSL_SUCCESS;
  430.     }
  431. }
  432.  
  433.  
  434. /* Update a QR factorisation for A= Q R ,  A' = A + u v^T,
  435.  
  436.  * Q' R' = QR + u v^T
  437.  *       = Q (R + Q^T u v^T)
  438.  *       = Q (R + w v^T)
  439.  *
  440.  * where w = Q^T u.
  441.  *
  442.  * Algorithm from Golub and Van Loan, "Matrix Computations", Section
  443.  * 12.5 (Updating Matrix Factorizations, Rank-One Changes)  
  444.  */
  445.  
  446. int
  447. gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
  448.               gsl_vector * w, const gsl_vector * v)
  449. {
  450.   const size_t M = R->size1;
  451.   const size_t N = R->size2;
  452.  
  453.   if (Q->size1 != M || Q->size2 != M)
  454.     {
  455.       GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
  456.     }
  457.   else if (w->size != M)
  458.     {
  459.       GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
  460.     }
  461.   else if (v->size != N)
  462.     {
  463.       GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
  464.     }
  465.   else
  466.     {
  467.       size_t j, k;
  468.       double w0;
  469.  
  470.       /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
  471.  
  472.          J_1^T .... J_(n-1)^T w = +/- |w| e_1
  473.  
  474.          simultaneously applied to R,  H = J_1^T ... J^T_(n-1) R
  475.          so that H is upper Hessenberg.  (12.5.2) */
  476.  
  477.       for (k = M - 1; k > 0; k--)
  478.     {
  479.       double c, s;
  480.       double wk = gsl_vector_get (w, k);
  481.       double wkm1 = gsl_vector_get (w, k - 1);
  482.  
  483.       create_givens (wkm1, wk, &c, &s);
  484.       apply_givens_vec (w, k - 1, k, c, s);
  485.       apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
  486.     }
  487.  
  488.       w0 = gsl_vector_get (w, 0);
  489.  
  490.       /* Add in w v^T  (Equation 12.5.3) */
  491.  
  492.       for (j = 0; j < N; j++)
  493.     {
  494.       double r0j = gsl_matrix_get (R, 0, j);
  495.       double vj = gsl_vector_get (v, j);
  496.       gsl_matrix_set (R, 0, j, r0j + w0 * vj);
  497.     }
  498.  
  499.       /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
  500.          Equation 12.5.4 */
  501.  
  502.       for (k = 1; k < GSL_MIN(M,N+1); k++)
  503.     {
  504.       double c, s;
  505.       double diag = gsl_matrix_get (R, k - 1, k - 1);
  506.       double offdiag = gsl_matrix_get (R, k, k - 1);
  507.  
  508.       create_givens (diag, offdiag, &c, &s);
  509.       apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
  510.  
  511.       gsl_matrix_set (R, k, k - 1, 0.0);    /* exact zero of G^T */
  512.     }
  513.  
  514.       return GSL_SUCCESS;
  515.     }
  516. }
  517.  
  518. int
  519. gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
  520. {
  521.   const size_t M = R->size1;
  522.   const size_t N = R->size2;
  523.  
  524.   if (M != N)
  525.     {
  526.       return GSL_ENOTSQR;
  527.     }
  528.   else if (Q->size1 != M || b->size != M || x->size != M)
  529.     {
  530.       return GSL_EBADLEN;
  531.     }
  532.   else
  533.     {
  534.       /* compute sol = Q^T b */
  535.  
  536.       gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
  537.  
  538.       /* Solve R x = sol, storing x in-place */
  539.  
  540.       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
  541.  
  542.       return GSL_SUCCESS;
  543.     }
  544. }
  545.